# ******** Generate Figures for the Paper *******
# Set up plot parameters
par_mar_setting = c(2.8, 2.8, 1.5, 1)  
mpg_setting  = c(1.5, 0.5, 0)
# Controls whether all all of the plots should output or not
output_to_files = TRUE  
# FigurePath is the path for figure output. The master file "Main_Step2.R" sets paths.

# --------- Figure 3a: Plot the Simulated Path of Bayesian and Diagnostic Beliefs  -------------
# parameter inputs -- use the diagnostic belief parameter set
lambda_H = 0.638
lambda_L = 0.001
lambda_H_to_L = 0.48
lambda_L_to_H = 0.11
# Then load the relevant functions for belief process. 
source("Belief_Process.R")
## Check the basic behavior of lambda
dt=1/12
result = as.data.table( simulation_hat_lambda(N_years = 100) )
result[ , lambda_theta1 := irrational_lambda(lambda, dt,  1, theta=1 )  ]
print("--------- Generating Figure 3(a) ----------")
pdf( paste0(FigurePath, "Belief Process/Belief_Process.pdf"), width = 5, height = 5)
par(mar = c(3.5, 3.5, 0.5, 0.5), mfrow=c(1,1))
MyPlot( yseries_input =  result[ i>38*12&i<50*12 , list((i-458)*dt, lambda, lambda_theta1)], LegendNames = c( "theta=0 (Bayesian Belief)", "theta=1" ), ylab = "perceived frequency of crises", xlab="years", ylim = c(lambda_L*0.9, lambda_H*1.2) )
abline( h=c(lambda_L, lambda_H), lty=3 )
dev.off()

# --------- Load the Simulations of Three Models -------------
option="main_paper"
print(">>>>>>>> Loading Simulated Data (This will take about 3 minutes) ")
source("Codes_Model/preprocessing_simulation_data.R")

#------------------- Generate the fitted data needed for Figure 3 (b) ------------------
unique_index = !duplicated( round(TB_Diagnostic$lambda_rational, digits=5) )
library(ConSpline)
TB_output = TB_Diagnostic[unique_index ,list(lambda_rational,  lambda_diagnostic=lambda )]
TB_output = TB_output[order(lambda_rational)]
write_xlsx(TB_output , "simulated data/lambda_relationship_diagnostic.xlsx" )


#---------------- Figure 6: Average Dynamcis Across Models ----------------------
# Average path around financial crises
gen_crises_dynamics <- function(TB, output_file_name=""){
  # Average values. 
  TB_dynamics_around_crises = TB[ !is.na(Years_around_crisis) , lapply(.SD, mean,na.rm=TRUE), by=Years_around_crisis ]
  TB_dynamics_around_crises = TB_dynamics_around_crises[ order(Years_around_crisis) ]
  TB_dynamics_around_crises[ , credit_spread_change := (credit_spread-mean(TB$credit_spread))/mean(TB$credit_spread) ]
  TB_dynamics_around_crises[ , bank_credit_change := (bank_credit-mean(bank_credit))/mean(bank_credit) ]
  TB_dynamics_around_crises[ , GDP_change := (productivity-mean(TB$productivity))/mean(TB$productivity) ]

  if( output_file_name!="" &&output_to_files){
    pdf(paste0(FigurePath,"dynamics_around_crises/",output_file_name), width = 6, height = 5)
    par(mar = c(3, 3, 1, 1), mfrow=c(2,2))
    MyPlot( TB_dynamics_around_crises$Years_around_crisis, TB_dynamics_around_crises$credit_spread_change , mgp = c(1.5, 0.5, 0),xlab = "years around crises", ylab="Credit Spread"); abline( h=0, lty=2 )
    MyPlot( TB_dynamics_around_crises$Years_around_crisis, TB_dynamics_around_crises$bank_credit_change , mgp = c(1.5, 0.5, 0),xlab = "years around crises", ylab="Bank Credit"); abline( h=0, lty=2 )
    MyPlot( TB_dynamics_around_crises$Years_around_crisis, TB_dynamics_around_crises$GDP_change, mgp = c(1.5, 0.5, 0),xlab = "years around crises", ylab="GDP" ); abline( h=0, lty=2 )
    dev.off()
  }
  par(mar = c(3, 3, 1, 1), mfrow=c(2,2))
  MyPlot( TB_dynamics_around_crises$Years_around_crisis, TB_dynamics_around_crises$credit_spread_change , mgp = c(1.5, 0.5, 0),xlab = "years around crises", ylab="Credit Spread"); abline( h=0, lty=2 )
  MyPlot( TB_dynamics_around_crises$Years_around_crisis, TB_dynamics_around_crises$bank_credit_change , mgp = c(1.5, 0.5, 0),xlab = "years around crises", ylab="Bank Credit"); abline( h=0, lty=2 )
  MyPlot( TB_dynamics_around_crises$Years_around_crisis, TB_dynamics_around_crises$GDP_change, mgp = c(1.5, 0.5, 0),xlab = "years around crises", ylab="GDP" ); abline( h=0, lty=2 )
  return(TB_dynamics_around_crises)
}
# For Slides
gen_crises_dynamics_slides <- function(TB_dynamics_around_crises, output_file_name=""){
  pdf(paste0(FigurePath,"dynamics_around_crises/",output_file_name), width = 7.5, height = 2.5)
  par(mar = c(3, 3, 1, 1), mfrow=c(1,3))
  MyPlot( TB_dynamics_around_crises$Years_around_crisis, TB_dynamics_around_crises$credit_spread_change , mgp = c(1.5, 0.5, 0),xlab = "years around crises", ylab="Credit Spread"); abline( h=0, lty=2 )
  MyPlot( TB_dynamics_around_crises$Years_around_crisis, TB_dynamics_around_crises$bank_credit_change , mgp = c(1.5, 0.5, 0),xlab = "years around crises", ylab="Bank Credit"); abline( h=0, lty=2 )
  MyPlot( TB_dynamics_around_crises$Years_around_crisis, TB_dynamics_around_crises$GDP_change, mgp = c(1.5, 0.5, 0),xlab = "years around crises", ylab="GDP" ); abline( h=0, lty=2 )
  dev.off()
}
TB_dynamics_Benchmark = gen_crises_dynamics( TB_benchmark )
TB_dynamics_Bayesian = gen_crises_dynamics(TB_Bayesian )
TB_dynamics_Diagnostic = gen_crises_dynamics(TB_Diagnostic)
if(output_to_files==TRUE){
  pdf(paste0(FigurePath,"dynamics_around_crises/all.pdf"), width = 9, height = 3)
}
par(mar = c(3, 3, 1, 1), mfrow=c(1,3))
MyPlot( TB_dynamics_Bayesian$Years_around_crisis, TB_dynamics_Bayesian$credit_spread_change , mgp = c(1.5, 0.5, 0),xlab = "years around crises", ylab="Credit Spread", ylim=c(-0.25, 0.6));  
abline( h=0, lty=2 )
lines( TB_dynamics_Diagnostic$Years_around_crisis,   TB_dynamics_Diagnostic$credit_spread_change, col="red", lty=2, lwd=2 )
lines( TB_dynamics_Benchmark$Years_around_crisis,   TB_dynamics_Benchmark$credit_spread_change, col="black", lty=2, lwd=1.5 )
legend( "topright", legend=c( "Bayesian Belief", "Diagnostic Belief", "Static Belief" ), bty="n", lty=c(1,2,3), col=c("blue","red","black"), lwd=rep(2,3) )

MyPlot( TB_dynamics_Bayesian$Years_around_crisis, TB_dynamics_Bayesian$bank_credit_change , mgp = c(1.5, 0.5, 0), xlab = "years around crises", ylab="Bank Credit", ylim=c(-1, 0.7)); 
abline( h=0, lty=2 )
lines( TB_dynamics_Diagnostic$Years_around_crisis,   TB_dynamics_Diagnostic$bank_credit_change, col="red", lty=2, lwd=2 )
lines( TB_dynamics_Benchmark$Years_around_crisis,   TB_dynamics_Benchmark$bank_credit_change, col="black", lty=2, lwd=1.5 )
legend( "topright", legend=c( "Bayesian Belief", "Diagnostic Belief", "Static Belief" ), bty="n", lty=c(1,2,3), col=c("blue","red","black"), lwd=rep(2,3) )

MyPlot( TB_dynamics_Bayesian$Years_around_crisis, TB_dynamics_Bayesian$GDP_change , mgp = c(1.5, 0.5, 0), xlab = "years around crises", ylab="GDP", ylim=c(-0.13, 0.1)); 
abline( h=0, lty=2 )
lines( TB_dynamics_Diagnostic$Years_around_crisis,   TB_dynamics_Diagnostic$GDP_change, col="red", lty=2, lwd=2 )
lines( TB_dynamics_Benchmark$Years_around_crisis,   TB_dynamics_Benchmark$GDP_change, col="black", lty=2, lwd=1.5 )
legend( "topright", legend=c( "Bayesian Belief", "Diagnostic Belief", "Static Belief" ), bty="n", lty=c(1,2,3), col=c("blue","red","black"), lwd=rep(2,3) )
if(output_to_files==TRUE){
  dev.off()
}


# --------------------- Figure 7: Volatility and Sharpe Ratio Dynamics around Crises ---------------------
pdf(paste0(FigurePath,"dynamics_around_crises/total_capital_return_empirical_vol.pdf"), width = 9, height = 3)
par(mar = c(3, 3, 1, 1), mfrow=c(1,3))
MyPlot( TB_dynamics_Benchmark$Years_around_crisis, TB_dynamics_Benchmark$capital_return_vol_simulation, mgp = c(1.5, 0.5, 0),xlab = "years around crises", ylab="vol(capital return) in simulation", main="Static Belief"); abline( h=0, lty=2 )
MyPlot( TB_dynamics_Bayesian$Years_around_crisis, TB_dynamics_Bayesian$capital_return_vol_simulation, mgp = c(1.5, 0.5, 0),xlab = "years around crises", ylab="vol(capital return) in simulation", main="Bayesian Belief"); abline( h=0, lty=2 )
MyPlot( TB_dynamics_Diagnostic$Years_around_crisis, TB_dynamics_Diagnostic$capital_return_vol_simulation, mgp = c(1.5, 0.5, 0),xlab = "years around crises", ylab="vol(capital return) in simulation", main="Diagnostic Belief"); abline( h=0, lty=2 )
dev.off()

pdf(paste0(FigurePath,"dynamics_around_crises/total_sharpe_ratio.pdf"), width = 9, height = 3)
par(mar = c(3, 3, 1, 1), mfrow=c(1,3))
MyPlot( TB_dynamics_Benchmark$Years_around_crisis, TB_dynamics_Benchmark$total_Sharpe, mgp = c(1.5, 0.5, 0),xlab = "years around crises", ylab="total Sharpe ratio", main="Static Belief"); abline( h=0, lty=2 )
MyPlot( TB_dynamics_Bayesian$Years_around_crisis, TB_dynamics_Bayesian$total_Sharpe, mgp = c(1.5, 0.5, 0),xlab = "years around crises", ylab="total Sharpe ratio", main="Bayesian Belief"); abline( h=0, lty=2 )
MyPlot( TB_dynamics_Diagnostic$Years_around_crisis, TB_dynamics_Diagnostic$total_Sharpe, mgp = c(1.5, 0.5, 0),xlab = "years around crises", ylab="total Sharpe ratio", main="Diagnostic Belief"); abline( h=0, lty=2 )
dev.off()
# Average total Sharpe ratio. 

pdf(paste0(FigurePath,"dynamics_around_crises/total_vol.pdf"), width = 9, height = 3)
par(mar = c(3, 3, 1, 1), mfrow=c(1,3))
MyPlot( TB_dynamics_Benchmark$Years_around_crisis, sqrt(TB_dynamics_Benchmark$total_risk_unadjusted), mgp = c(1.5, 0.5, 0),xlab = "years around crises", ylab="total risks", main="Static Belief"); abline( h=0, lty=2 )
abline(h=1, lty=3)
MyPlot( TB_dynamics_Bayesian$Years_around_crisis, sqrt(TB_dynamics_Bayesian$total_risk_unadjusted), mgp = c(1.5, 0.5, 0),xlab = "years around crises", ylab="total risks", main="Bayesian Belief"); abline( h=0, lty=2 )
abline(h=1, lty=3)
MyPlot( TB_dynamics_Diagnostic$Years_around_crisis, sqrt(TB_dynamics_Diagnostic$total_risk_unadjusted), mgp = c(1.5, 0.5, 0),xlab = "years around crises", ylab="total risks", main="Diagnostic Belief"); abline( h=0, lty=2 )
abline(h=1, lty=3)
dev.off()




# ------------------------ Figure 8: Stationary Distribution of Stata Variables in the Model ---------------------------
plot_density <- function(x, adjust=0, xlab="x", bw=0.01, from=0, to=1, file_name){
  if(output_to_files){
    pdf( paste0("../Figures/stationary distribution/", file_name), width = 5, height = 4)
  }
  par(mar = c(3,3,1.5,1) , mfrow=c(1,1))
  plot(  density( x+adjust,  bw=bw, from=from, to=to ), col="blue", lwd=2, main="", xlab=xlab,   ylab="density", mgp=mpg_setting  )
  if(output_to_files){
    dev.off()
  }
  par(mar = c(3,3,1.5,1) , mfrow=c(1,1))
  plot(  density( x+adjust,  bw=bw, from=from, to=to ), col="blue",lwd=2,main="", xlab=xlab,   ylab="density", mgp=mpg_setting  )
}
plot_density(  TB_benchmark$w, adjust=0.01,  xlab="w",  bw=0.01,  from=0,  to=0.4, file_name = "density_benchmark_w.pdf" )
plot_density(  TB_Bayesian$w, adjust=0.01,  xlab="w",  bw=0.01,  from=0,  to=0.4, file_name = "density_Bayesian_w.pdf" )
plot_density(  TB_Diagnostic$w, adjust=0.03,  xlab="w",  bw=0.01,  from=0,  to=0.4, file_name = "density_Diagnostic_w.pdf" )
plot_density(  TB_benchmark$lambda, adjust=0,  xlab="lambda",  bw=0.00001,  from=mean(TB_benchmark$lambda)*0.98,  to=mean(TB_benchmark$lambda)*1.02, file_name = "density_benchmark_lambda.pdf" )
plot_density(  TB_Bayesian$lambda, adjust=0,  xlab="lambda",  bw=0.01,  from=0,  to=0.5, file_name = "density_Bayesian_lambda.pdf" )
plot_density(  TB_Diagnostic$lambda, adjust=0.01,  xlab="lambda",  bw=0.01,  from=0,  to=0.72, file_name = "density_Diagnostic_lambda.pdf" )
plot_density_all <- function( TB_states, legend_position="topleft",  xlab="x", bw=0.01, from=0, to=0.4, ylim =c(0,5) , file_name=""){
  if(output_to_files && file_name!=""){
    pdf( paste0("../Figures/stationary distribution/", file_name), width = 5, height = 4.5)
  }
  par(mar = c(3,3,1.5,1) , mfrow=c(1,1))
  plot(  density( TB_states[,1],  bw=bw, from=from, to=to ), col="blue", lwd=2, main="", xlab=xlab,   ylab="density", mgp=mpg_setting, ylim= ylim)
  lines( density( TB_states[,2]+0.03,  bw=bw, from=from, to=to ), col="red", lwd=2, lty=2  )
  lines( density( TB_states[,3],  bw=bw, from=from, to=to ), col="purple", lwd=2, lty=2  )
  legend( legend_position, legend = c( "Bayesian Belief", "Diagnostic Belief",  "Static Belief"), bty = "n", lty= c(1,2,3) , col=c("blue","red","purple"),  lwd=c(2,2,2)  )
  if(output_to_files && file_name!=""){
    dev.off()
  }
}
TB_states = cbind( TB_Bayesian$w, TB_Diagnostic$w, TB_benchmark$w )
plot_density_all(  TB_states,  xlab="w",  bw=0.01,  from=0,  to=0.35, ylim=c(0,11) )
TB_states = cbind( TB_Bayesian$lambda, TB_Diagnostic$lambda, TB_benchmark$lambda-0.02 )
plot_density_all(  TB_states, legend_position="topright",  xlab="lambda",  bw=0.02,  from=0,  to=0.75, ylim=c(0,20))


# ------------------ Figure 9: Histogram of GDP growth for the next three years ------------------
TB_Data = read_xlsx(paste0(FigurePath,"figures_from_data_of_Krishnamurthy_and_Muir/Histogram data.xlsx"),range="B1:B43")
names(TB_Data)<-"GDP_drop"
h <- hist(TB_Data$GDP_drop/100, breaks = 20, plot=FALSE) 
multiplier=1.1
DT_density = density( TB_Bayesian[crisis==1]$GDP_3Y_change )
h$counts=h$counts/max(h$counts)*max(DT_density$y)*multiplier
# Then plot the data 
plot_GDP_growth_density <- function(file_name, GDP_drop, model_name, bw=0.02, xlab="3 Year GDP Growth"){  
  for(i in 1:2){
    if(i==2&&output_to_files&&file_name!=""){
      pdf( file_name, width=7, height=4 )
    }
    DT_density = density( TB[crisis==1]$GDP_3Y_change )
    par(mar = c(3, 3, 1, 1), mfrow=c(1,1))
    h <- hist(TB_Data$GDP_drop/100, breaks = 25, plot=FALSE) 
    multiplier=1.3
    h$counts=h$counts/max(h$counts)*max(DT_density$y)*multiplier
    plot(h , main="", xlab=xlab, ylab="Density", mgp = c(2, 0.8, 0), col="lightgoldenrod3",  xlim=c(-0.4, 0.1), ylim=c(0, max(DT_density$y)*multiplier) )
    lines( density(GDP_drop, bw=bw ), col="black", lwd=2, lty=1 )
    legend(  "topleft", legend= c("Data", paste0(model_name," Model") ) , bty="n", col=c("lightgoldenrod3", "black"), lwd = c(14,2)   )
    if(i==2&&output_to_files&&file_name!=""){
      dev.off()
    }
  }
}
if(output_to_files){
  pdf( "../Figures/gdp_growth_around_crises/distributions.pdf" , width=8, height=5.5 )
}
par(mar = c(3, 3, 1, 1), mfrow=c(1,1))
plot(h , main="", xlab="3 Year GDP Growth", ylab="Density", mgp = c(2, 0.8, 0), col="lightgoldenrod3",  xlim=c(-0.4, 0.1), ylim=c(0, max(DT_density$y)*multiplier) )
lines( density(TB_Bayesian[crisis==1]$GDP_3Y_change, bw=0.02 ), col="blue", lwd=2, lty=1 )
lines( density(TB_Diagnostic[crisis==1]$GDP_3Y_change, bw=0.02 ), col="red", lwd=2, lty=2 )
lines( density(TB_benchmark[crisis==1]$GDP_3Y_change, bw=0.02 ), col="black", lwd=2, lty=3 )
legend(  "topleft", legend= c("Data", "Bayesian Belief", "Diagnostic Belief", "Static Belief" ) , bty="n", col=c("lightgoldenrod3", "blue", "red", "black"), lwd = c(14,2,2,2), lty=c(1,1,2,3)   )
if(output_to_files){
  dev.off()
}



# -----------  Preparations to Plot Conditional GDP Growth ---------------------
cut_off = c(0.3, 0.9)
cut_off_reverse = 1-cut_off[c(2,1)]
plot_conditional_density <- function( quantile_lower=0, quantile_upper = 0.2, conditioning_var_name = "credit_spread", y_var_name = "GDP_yearly_change",
                                      file_name = "../Figures/gdp_growth_around_crises/GDP_growth_density_on_low_credit_spread.pdf",
                                     xlab = "Three-Year GDP Growth", bw=0.008, xlim=NULL, ylim_upper=NA, plot_option=1, plot_vertical_line=TRUE,
                                     width=7, height=4){
  y_quantity1 = TB_benchmark[[y_var_name]]
  y_quantity2 = TB_Bayesian[[y_var_name]]
  y_quantity3 = TB_Diagnostic[[y_var_name]]
  conditional_quantity1 = TB_benchmark[[conditioning_var_name]]
  conditional_quantity2 = TB_Bayesian[[conditioning_var_name]]
  conditional_quantity3 = TB_Diagnostic[[conditioning_var_name]]
  for(i in 1:2){
    if(i==2&&output_to_files){
      pdf( file_name , width=width, height=height )
    }
    par(mar = c(3, 3, 1, 1), mfrow=c(1,1))
    lty_vec = c(2,1,4);  col_vec = c("blue","red","seagreen");   lwd_vec = c(2,1.5,2.5)
    legend_names = c("Static Belief","Bayesian","Diagnostic")
    y1 = y_quantity1[conditional_quantity1<=quantile(conditional_quantity1,quantile_upper)&conditional_quantity1>=quantile(conditional_quantity1,quantile_lower)]
    y2 = y_quantity2[conditional_quantity2<=quantile(conditional_quantity2,quantile_upper)&conditional_quantity2>=quantile(conditional_quantity2,quantile_lower)]
    y3 = y_quantity3[conditional_quantity3<=quantile(conditional_quantity3,quantile_upper)&conditional_quantity3>=quantile(conditional_quantity3,quantile_lower)]
    density1=density(y1,  bw=bw);  density2=density(y2,  bw=bw);  density3=density(y3,  bw=bw);
    ylim = c( 0,  max( density1$y, density2$y, density3$y, ylim_upper, na.rm = TRUE )  )
    if(is.null(xlim)){
      plot(  density2,   lwd=lwd_vec[2], col=col_vec[2], lty=lty_vec[2], mgp = c(2, 0.8, 0), main="", xlab=xlab, ylim=ylim )
    }else{
      plot(  density2,   lwd=lwd_vec[2], col=col_vec[2], lty=lty_vec[2], mgp = c(2, 0.8, 0), main="", xlab=xlab, ylim=ylim, xlim=xlim )
    }
    if(plot_option<=2){
      if(plot_option==1){
        index = c(2,1)
        y_quantity = c(y_quantity2,y_quantity1)
        lines( density1 ,  lwd=lwd_vec[1], col=col_vec[1], lty=lty_vec[1])
      }else if(plot_option==2){
        index = c(2,3)
        y_quantity = c(y_quantity2,y_quantity3)
        lines( density3 ,  lwd=lwd_vec[3], col=col_vec[3], lty=lty_vec[3])
      }
      if(plot_vertical_line){
        abline( v=quantile(c(y_quantity1,y_quantity2,y_quantity3), 0.04), lty=3, lwd=1.5 )
      }
      legend(  "topleft", legend=legend_names[index], bty="n", lty=lty_vec[index], col=col_vec[index], lwd = lwd_vec[index]  )
    }else{  # Special treatment for all plots
      index = c(2,3,1)
      plot(  density2,   lwd=2, col="blue", lty=1, mgp = c(2, 0.8, 0), main="", xlab=xlab, ylim=ylim, xlim=xlim )
      y_quantity = c(y_quantity1,y_quantity2,y_quantity3)
      lines( density3 ,  lwd=2, col="red", lty=2)
      lines( density1 ,  lwd=2.2, col="purple", lty=3)
      if(plot_vertical_line){
        abline( v=quantile(c(y_quantity1,y_quantity2,y_quantity3), 0.04), lty=3, lwd=1.5 )
      }
      legend(  "topleft", legend=legend_names[index], bty="n", lty=c(1,2,3), col=c("blue","red","purple"), lwd = c(2,2,2.2)  )
    }
    if(i==2&&output_to_files){
      dev.off()
    }
  }
}
#------------------------- Figure 11: Density of Next-Year GDP Growth Conditional on Bank Credit/GDP -----------------------
y_var_name = "GDP_yearly_change";  xlab="Next-Year GDP Growth";
plot_conditional_density(quantile_lower=0, quantile_upper = cut_off_reverse[1], conditioning_var_name = "bank_credit_GDP", y_var_name = y_var_name,
                         file_name = "../Figures/gdp_growth_around_crises/GDP_growth_density_on_low_bank_credit.pdf",
                         xlab = xlab, xlim=c(-0.25,0.15), ylim_upper = 18.5, bw=0.008)
plot_conditional_density(quantile_lower=cut_off_reverse[2], quantile_upper = 1, conditioning_var_name = "bank_credit_GDP", y_var_name = y_var_name,
                         file_name = "../Figures/gdp_growth_around_crises/GDP_growth_density_on_high_bank_credit.pdf",
                         xlab = xlab, xlim=c(-0.25,0.15), ylim_upper = 18.5, bw=0.008)
plot_conditional_density(quantile_lower=0, quantile_upper = cut_off_reverse[1], conditioning_var_name = "bank_credit_GDP", y_var_name = y_var_name,
                         file_name = "../Figures/gdp_growth_around_crises/GDP_growth_density_on_low_bank_credit_Diagnostic.pdf",
                         xlab = xlab, xlim=c(-0.25,0.2), ylim_upper = 18.5, bw=0.012, plot_option = 2)
plot_conditional_density(quantile_lower=cut_off_reverse[2], quantile_upper = 1, conditioning_var_name = "bank_credit_GDP", y_var_name = y_var_name,
                         file_name = "../Figures/gdp_growth_around_crises/GDP_growth_density_on_high_bank_credit_Diagnostic.pdf",
                         xlab = xlab, xlim=c(-0.25,0.2), ylim_upper = 18.5, bw=0.008, plot_option = 2)

#------------------------- Figure 12: Density of Next-Year GDP Growth Conditional on Credit Spread ------------------
y_var_name = "bank_equity_excess_return"; xlab="One-Year Bank Equity Return";
plot_conditional_density(quantile_lower=0, quantile_upper = cut_off[1], conditioning_var_name = "bank_credit_GDP", y_var_name = y_var_name,
                         file_name = "../Figures/banker_wealth_growth/w_growth_density_on_low_bank_credit.pdf",
                         xlab = xlab, xlim=c(-0.8,0.7), ylim_upper = 7, bw=0.04, plot_vertical_line = FALSE)
plot_conditional_density(quantile_lower=cut_off[2], quantile_upper = 1, conditioning_var_name = "bank_credit_GDP", y_var_name = y_var_name,
                         file_name = "../Figures/banker_wealth_growth/w_growth_density_on_high_bank_credit.pdf",
                         xlab = xlab, xlim=c(-0.8,0.7), ylim_upper = 7, bw=0.03, plot_vertical_line = FALSE)

#------------------------- Figure 15: Density of Next-Year GDP Growth Conditional on Bank Credit/GDP -----------------------
y_var_name = "GDP_yearly_change";  xlab="Next-Year GDP Growth";
plot_conditional_density(quantile_lower=0, quantile_upper = cut_off[1], conditioning_var_name = "credit_spread", y_var_name = y_var_name,
                         file_name = "../Figures/gdp_growth_around_crises/GDP_growth_density_on_low_credit_spread.pdf",
                         xlab = xlab, xlim=c(-0.25,0.2), ylim_upper = 18.5, bw=0.008)
plot_conditional_density(quantile_lower=cut_off[2], quantile_upper = 1, conditioning_var_name = "credit_spread", y_var_name = y_var_name,
                         file_name = "../Figures/gdp_growth_around_crises/GDP_growth_density_on_high_credit_spread.pdf",
                         xlab = xlab, xlim=c(-0.25,0.2), ylim_upper = 18.5, bw=0.02)
plot_conditional_density(quantile_lower=0, quantile_upper = cut_off[1], conditioning_var_name = "credit_spread", y_var_name = y_var_name,
                         file_name = "../Figures/gdp_growth_around_crises/GDP_growth_density_on_low_credit_spread_Diagnostic.pdf",
                         xlab = xlab, xlim=c(-0.25,0.2), ylim_upper = 18.5, bw=0.008, plot_option = 2)
plot_conditional_density(quantile_lower=cut_off[2], quantile_upper = 1, conditioning_var_name = "credit_spread", y_var_name = y_var_name,
                         file_name = "../Figures/gdp_growth_around_crises/GDP_growth_density_on_high_credit_spread_Diagnostic.pdf",
                         xlab = xlab, xlim=c(-0.25,0.2), ylim_upper = 18.5, bw=0.02, plot_option = 2)

